**********************************************
*** 	       Time Series Plots          ****
**********************************************
global TAXVARS uscetr_soi paystax
global SUBSAMPLE firmtype_alt salq tincq
	
foreach insample in $SUBSAMPLE {
	foreach v in $TAXVARS  {
		
		*Hack so Stata doesn't "forget" the v local after the preserve/restore
		global v `v'
		
		*Read data
		u "$cur_sample.dta", clear
		do $cleaning/renameSOI
		
		g pos_inc = pidom_soi>0
		
		*Winsorize to improve data visualizations
		winsor2 uscetr_soi, cuts(0 95) replace //by(year)
		replace uscetr_soi=100*uscetr_soi
		*Yearly Mean
		gegen `v'_overall = mean(`v') [pw=weightSOI], by(year)
					
		*Total taxes, employment, firms in each group
		g xx = 1
		gegen tax_overall = sum(tax) [pw=weightSOI], by(year)
		gegen taxx = sum(tax) [pw=weightSOI], by(year `insample')
		replace tax = taxx
		drop taxx
		gegen emp_overall = sum(emp_w2) [pw=weightSOI], by(year)
		gegen emp = sum(emp_w2) [pw=weightSOI], by(year `insample')
		gegen sal_overall = sum(revenue) [pw=weightSOI], by(year)
		gegen sal = sum(revenue) [pw=weightSOI], by(year `insample')		
		
		gegen N_overall = sum(xx) [pw=weightSOI], by(year)
		gegen N = sum(xx) [pw=weightSOI], by(year `insample')
		drop xx
		
		*For ETRS, drop firms for which we can't compute it
		if inlist("`v'","uscetr_soi") {
			drop if pidom_soi<=0
		}
		
		*Total number of firms 
		gdistinct firmid
		local N = string(r(ndistinct),"%9.0fc")
		
		*Number of firms in each subgroup
		forvalues f = 1/4{
			gdistinct firmid if firmtype_alt==`f'
			local N_`f' = string(r(ndistinct),"%9.0fc")
		}	
		
		*ETR
		g etr = uscetr_soi
		g pays = paystax*100
		g etr_alt = uscetr_soi_alt
		g pos = pos_inc*100

		foreach j in etr pays etr_alt pos{
		foreach i in 2016 2019 {	
			sum `j' if year==`i' [aw=weightSOI]
			local `j'_`i' = string(r(mean),"%9.1f")
		}
		}		
		
		*Averages by year and group
		gcollapse (mean) *_overall `v' etr pays etr_alt pos tax emp sal N [pw=weightSOI], ///
				by(year `insample') labelformat(#sourcelabel#) fast
				
		*Reshape for line charts		
		local v $v
		reshape wide `v' tax emp N sal, i(year *_overall) j(`insample')
		
		*Shares in total # firms; employment; federal taxes
		g taxw = .
		g empw = . 
		g Nw = .
		g salw = .
		foreach j in tax emp N sal {
			forval f = 1/5 {
				cap replace `j'w = `j'`f'/`j'_overall if mi(`j'w)
			}
		}
		
		
		**Plot title,Y-axis label and format
		local minyr 2010
		local endyr 2020
		keep if year>=`minyr' & year<=`endyr'
	
		if strpos("`v'","cetr")>0 {
			local format "%9.0f"
			local title "US Cash ETR"
			local xlab "Percent"
			local fmt "%9.0f"
		}
		if strpos("`v'","paystax")>0 {
			local format "%9.0f"
			local title "Share of Firms Paying >0 Federal Income Tax"
			local xlab "Share"
			local fmt "%9.2f"
		}
		if strpos("`v'","pos_inc")>0 {
			local format "%9.0f"
			local title "Share of Firms with Positive Book Profit"
			local xlab "Share"
			local fmt "%9.2f"
		}
		
		if inlist("`insample'","firmtype_alt") {
			local lab1 "Pub MNE"
			local lab2 "Pub Dom"
			local lab3 "Priv MNE"
			local lab4 "Large Priv Dom"
			local lab5 "Small Priv Dom"


		}
		if inlist("`insample'","tincq") {
			local title `"`title'" "by Centile of Firm Taxable Income Distribution"'
			local lab1 "Bottom 90%"
			local lab2 "90-95%"
			local lab3 "95%+"


		}
		if inlist("`insample'","salq") {
			local title `"`title'" "by Centile of Firm Gross Receipts Distribution"'
			local lab1 "Bottom 90%"
			local lab2 "90-95%"
			local lab3 "95%+"


		}
		
	foreach yt in 2019 2020 {

		if inlist("`insample'", "firmtype_alt") {
			*Plot
			twoway  (scatter `v'_overall `v'1 `v'2 `v'3 `v'4 `v'5 year if year<=`yt', connect(l l l l l l) ///
				lpattern(solid dash solid shortdash solid longdash) msymbol(none none diamond none x none)	///
				lwidth(thick medium medium thick medium medium) msize(small small medium small huge small) ///
				legend(order(1 "Overall" 2 "`lab1'" 3 "`lab2'" 4 "`lab3'" 5 "`lab4'" 6 "`lab5'") row(2) pos(12) size(medsmall)) ///
				xaxis(1 2) xlab(`minyr'(1)`yt', axis(1) labsize(small)) xlab(`minyr' "`xlab'", noticks axis(2) labsize(medsmall)) ) ///
				(scatter `v'1 `v'2 `v'3 `v'4 `v'5 year if year==2016 [w=taxw], mcolor(purple purple purple purple purple purple) m(Oh Oh Oh Oh Oh 0h)), ///
				ytitle("") xtitle("") xtitle("", axis(2)) xsize(4.5) ysize(4) ///
				ylab(,format(`fmt') labsize(medsmall)) 				
		}		
		
		if inlist("`insample'", "salq","tincq"){
			*Plot
			twoway  (scatter `v'_overall `v'1 `v'2 `v'3 year if year<=`yt', ///
				connect(l l l l) lpattern(solid dash longdash solid)	///
				msymbol(none none none diamond) lwidth(thick medium medium medium) ///
				legend(order(1 "Overall" 2 "`lab1'" 3 "`lab2'" 4 "`lab3'" /*4 "`lab4'" 5 "`lab5'"*/) row(1) pos(12) size(medsmall)) ///
				xaxis(1 2) xlab(`minyr'(1)`yt', axis(1) labsize(small)) xlab(`minyr' "`xlab'", noticks axis(2) labsize(medsmall)) ) ///
				(scatter `v'1 `v'2 `v'3 year if year==2016 [w=taxw], mcolor(purple purple purple purple purple) m(Oh Oh Oh Oh Oh)), ///
				ytitle("") xtitle("") xtitle("", axis(2)) xsize(4.5) ysize(4) ///
				ylab(,format(`fmt') labsize(medsmall)) 
		}
		
		graph export "$results/`v'_`insample'_overall_`yt'.pdf", replace
		*graph export "$results/`v'_`insample'_overall_`yt'.png", replace

		}
	}
}
